#include "stdafx.h"

int			**grid;
double		boundaryHtc[NUM_MATERIAL][NUM_MATERIAL];
dieNum2D	dieNum[NUM_DIE];
plateNum2D	plateNum;
int			airThickNum;
int			row, col;
double		**tempCurr, **tempNext;

int initBoundaryHtc(material substance[], double dx)
{
	double r, k;
	for(int i = 0; i < NUM_MATERIAL; i++)
		for(int j = i; j < NUM_MATERIAL; j++)
		{
			if (substance[i].both[j] == 0) r = 0; else r = 1 / substance[i].both[j];
			if (i == 0 && (j == 1 || j == 2)) k = TOUCH_FACTOR; else k = 1;
			boundaryHtc[i][j] = boundaryHtc[j][i] = 
				k * dx / (r + dx / 2 / substance[i].htc + dx / 2 / substance[j].htc);
			if (i == 0 && (j == 1 || j == 2)) boundaryHtc[i][j] = boundaryHtc[j][i] = 7;
		}
	return 0;
}

dieNum2D getDieNum(die2D die, double dx)
{
	dieNum2D dieNum;
	dieNum.heightNum = (int)(die.height /dx);
	dieNum.widthNum = (int)(die.width / dx);
	return dieNum;
}

plateNum2D getPlateNum(plate2D plate, double dx)
{
	plateNum2D plateNum;
	plateNum.heightNum = (int)(plate.height /dx);
	plateNum.widthNum = (int)(plate.width / dx);
	return plateNum;
}

int initBlock(int up, int down, int left, int right, int id, double initTemp)
{
	for(int i = up; i < down; i++)
		for(int j = left; j < right; j++)
		{
			grid[i][j] = id;
			tempCurr[i][j] = tempNext[i][j] = initTemp;
		}
	return 0;
}

int initGrid(material substance[])
{
	grid = (int **)malloc(sizeof(int *) * row);
	tempCurr = (double **)malloc(sizeof(double *) * row);
	tempNext = (double **)malloc(sizeof(double *) * row);
	for(int i = 0; i < row; i++)
	{
		grid[i] = (int *)malloc(sizeof(int) * col);
		tempCurr[i] = (double *)malloc(sizeof(double) * col);
		tempNext[i] = (double *)malloc(sizeof(double) * col);
	}
	
	int up, down, left, right, k;
	//init block - Gas - Air
	k = 3;
	up = 0; down = row; left = 0; right = col; 
	initBlock(up, down, left, right, substance[k].id, substance[k].initTemp);

	//init block - Plate - 
	k = 0;
	up = airThickNum + dieNum[0].heightNum; down = up + plateNum.heightNum; 
	left = (col - plateNum.widthNum) / 2; right = left + plateNum.widthNum; 
	initBlock(up, down, left, right, substance[k].id, substance[k].initTemp);
	
	//init block - Die Up -  
	k = 1;
	up = airThickNum; down = up + dieNum[0].heightNum; 
	left = (col - dieNum[0].widthNum) / 2; right = left + dieNum[0].widthNum;
	initBlock(up, down, left, right, substance[k].id, substance[k].initTemp);
	
	//init block - Die Down - 
	k = 2;
	up = airThickNum + dieNum[0].heightNum + plateNum.heightNum; down = up + dieNum[1].heightNum; 
	left = (col -  dieNum[1].widthNum) / 2; right = left + dieNum[1].widthNum; 
	initBlock(up, down, left, right, substance[k].id, substance[k].initTemp);	

	tempCurr[0][0] = tempNext[0][0] = substance[0].initTemp;

	return 0;
}

double getMaxTimeStep(material substance[], double dx)
{
	double maxTime = 100000, temp;
	double dxdx = dx * dx;
	for(int i = 1; i < row - 1; i++)
		for(int j = 1; j < col - 1; j++)
		{
			temp = (substance[grid[i][j]].density * substance[grid[i][j]].sh * dxdx)
				/ (boundaryHtc[grid[i][j]][grid[i - 1][j]] + boundaryHtc[grid[i][j]][grid[i + 1][j]]
				+ boundaryHtc[grid[i][j]][grid[i][j - 1]] + boundaryHtc[grid[i][j]][grid[i][j + 1]]);
			if (temp < maxTime)
				maxTime = temp;
		}
	return maxTime;
}

int calNextTemp(material substance[], int i, int j, double timeStep, double dxdx)
{
	int k = grid[i][j];
	tempNext[i][j] = tempCurr[i][j] + 
		timeStep / substance[k].density / substance[k].sh / dxdx * 
		(boundaryHtc[k][grid[i][j - 1]]*(tempCurr[i][j - 1] - tempCurr[i][j])+
		boundaryHtc[k][grid[i][j + 1]]*(tempCurr[i][j + 1] - tempCurr[i][j])+
		boundaryHtc[k][grid[i - 1][j]]*(tempCurr[i - 1][j] - tempCurr[i][j])+
		boundaryHtc[k][grid[i + 1][j]]*(tempCurr[i + 1][j] - tempCurr[i][j]));
	return 0;
}

double tempFieldIteration(material substance[], double timeStep, double dx)
{
	double	dxdx = dx * dx;
	int		pointX = airThickNum + dieNum[0].heightNum + plateNum.heightNum / 2;
	int		pointY = col / 2;
	double	**tag;
	double	totalTime = 0;
	//output(tempCurr, row, col, totalTime, pointX, pointY);
	double	tempGoal = GOAL_TEMPERATURE;//substance[2].initTemp * 1.2;
	while(tempCurr[pointX][pointY] > tempGoal)
	{
		totalTime += timeStep;
		for(int i = 1; i < row - 1; i++)
			for(int j = 1; j < col -1; j++)
			{
				calNextTemp(substance, i, j, timeStep, dxdx);
				//tempCurr[i][j] = tempNext[i][j];
			}
		tag = tempCurr; tempCurr = tempNext; tempNext = tag;	
/*		if ((int)(totalTime / OUTPUT_FREQUENCY) > (int)((totalTime - timeStep) / OUTPUT_FREQUENCY)) 
			output(tempCurr, row, col, totalTime, pointX, pointY);*/		
	}
	//output(tempCurr, row, col, totalTime + OUTPUT_FREQUENCY * 2, pointX, pointY);
	return totalTime;
}

//********************************************************************************
double getCostTime(LARGE_INTEGER start, LARGE_INTEGER end)
{
	LARGE_INTEGER nFrequency;
	QueryPerformanceFrequency(&nFrequency);
	double time = (double)TIME_ACCURACY * (end.QuadPart - start.QuadPart) / nFrequency.QuadPart;
	return time / TIME_ACCURACY;
}

double fourier2D(struct material substance[], struct die2D die[], struct plate2D plate, double airThick, double dx)
{
	//timer
	LARGE_INTEGER	nStartCounter, nStopCounter;
	double			nTime;
	//timer starts
	QueryPerformanceCounter(&nStartCounter);

	initBoundaryHtc(substance, dx);

	dieNum[0] = getDieNum(die[0], dx);
	dieNum[1] = getDieNum(die[1], dx);
	plateNum = getPlateNum(plate, dx);
	airThickNum = (int)(airThick / dx);

	row = dieNum[0].heightNum + dieNum[1].heightNum + plateNum.heightNum + airThickNum * 2;
	col = (dieNum[0].widthNum > dieNum[1].widthNum ? dieNum[0].widthNum : dieNum[1].widthNum) + 2 * airThickNum;
	//printf("row = %d \t col = %d \n", row, col);

	initGrid(substance);

	//time step
	double maxTimeStep = getMaxTimeStep(substance, dx);
	//printf("Max Time Step = %f \n", maxTimeStep);
	double timeStep = TIME_STEP_FACTOR * maxTimeStep;

	double totalTime = tempFieldIteration(substance, timeStep, dx);

	//timer2 ends
	QueryPerformanceCounter(&nStopCounter);
	//get time
	nTime = getCostTime(nStartCounter, nStopCounter);

	printf("totalTime = %lf\n", totalTime);
	printf("totalStep = %d\n", (int)(totalTime / timeStep));
	printf("Compute Cycle = %f\n", nTime);

	return totalTime;
}